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ABSTRACT 

We examine the effects of hydrodynamics on the late stage kinetics in spinodal 
decomposition. From computer simulations of a lattice Boltzmann scheme we observe, 
for critical quenches, that single phase domains grow asymptotically like t", with a ~ 
.66 in two dimensions and a ~ 1.0 in three dimensions, both in excellent agreement 
with theoretical predictions. 
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The understanding of phase segregation kinetics has been enhanced over the last 
two decades by the development of dynamic renormalization group methods, an in- 
crease in computer resources, and a host of experimental results with which to com- 
pare 0. Despite this progress, there remain many open questions, especially 
regarding the growth of single phase domains and the scaling properties of the corre- 
lation or structure functions. Since phase segregation, either spinodal decomposition 
or nucleation, is a highly nonlinear process, theories rely mostly on approximation 
schemes, and exact results have been obtained only in certain limiting cases ^ . Sim- 
ulations designed to address these issues typically require extensive computational 
effort and have not, in our opinion, provided conclusive answers to many of the most 
fundamental questions. 

Phase segregating systems fall into two classes: those with hydrodynamic interac- 
tions (fluids) and those without (binary alloys, glasses). The latter class, has received 
far more attention both by computer simulation and by theory [||, 0. Only recently 
have there been attempts to carry out computer simulations of phase segregating sys- 
tems with hydrodynamic interactions. These include molecular dynamics (MD) simu- 
lations direct numerical simulation of time-dependent Ginzburg-Landau (TDGL) 
equations 0, cell-dynamical systems (CDS) 0, 0] and lattice gas (LG) H, Q and 
lattice Boltzmann (LB) models |T0[| . 

While MD simulations accurately represent the dynamics of real fluids, they are 
computationally demanding and at present may not be able to access the late stage 
(scaling) regime of spinodal decomposition. The TDGL approach eliminates the ex- 
act Newtonian particle dynamics in favor of a stochastic evolution governed by a 
coarse-grained phenomenological free energy functional. One then solves the con- 
comitant system of Langevin equations which couple order parameter fluctuations 
to hydrodynamic currents. This approach requires extensive numerical computation 
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and sacrifices both the Navier-Stokes behavior of the fluid and the interface phase 
dynamics. CDS methods further abstract the phase segregation process by replacing 
the Langevin equations with a much less numerically intensive dynamical map for the 
order parameter. The hydrodynamics they incorporate, however, have either been ap- 
proximated by an Oseen tensor or by a coupling of the velocity field to pressure and 
order parameter fluctuations. 

LG and LB models, on the other hand, provide alternative computational envi- 
ronments with which to study hydrodynamic phase segregation phenomena without 
the introduction of ad hoc relations between the order parameter fluctuations and 
the fluid dynamics. What makes these schemes so appealing is the natural way in 
which they can simulate the fluid properties, the phase segregation and the interface 
dynamics simultaneously. One drawback in the model presented here is that fluctu- 
ations only appear in the initial conditions. To make LB models more realistic tools 
with which to analyze phase transition kinetics probably involves the introduction of 
some type of thermal noise. 

Several groups ^ ^ have used LG and LB models to simulate phase segre- 
gation phenomena, but their work focused primarily on the qualitative features and 
did not address the details of dynamical scaling of the structure function |Tl|]. In 



this paper we focus on domain growth kinetics and scaling properties of the phase 
segregation process. We present the results of large-scale numerical simulations of 
one such immiscible LB model and make quantitative comparisons with theory. 

The LB method is a discrete, in space and time, microscopic, kinetic equation 
description for the evolution of the particle distribution function of a fluid |jl2[. The 



scheme described here is a modifled version of the immiscible fluid model introduced 



by Gunstensen et al [jlO| . In this model the fluid has two components represented, for 



example, by the colors red and blue. The microscopic dynamics of the particle distri- 



bution function consists of four steps: (1) free streaming, (2) collision, (3) interface 
perturbation and (4) recoloring. 

In free streaming the fluid (both red and blue components) moves to neighboring 
sites along the links of the underlying lattice. During the collision step these densities 
then relax toward a local equilibrium state. In LB schemes one is free to specify the 
local equilibrium state, and the particular choice for this state is one which leads to 
the Navier-Stokes equations in the long wavelength and low frequency limit. 

Let /j(x, t), //"(x, t) and /^^'(x, t) be the distribution functions for the total fluid, 
red (r) fluid and blue (b) fluid, respectively, at site x and time t moving along link in 
the i direction. Here /i(x, t) = fii^x, t) + (x, t), where i = 0, 1, ■ ■ ■, iV, and where N 
is the number of velocity states at each site. The 2 = state represents the portion 
of the fluid at rest. The kinetic or "lattice Boltzmann" equation for /j is written 

/,(x + e„ t + 1) - /,(x, t) = (]^(x, t) + (]f (x, t), (1) 

where VLl is the term representing the rate of change of fi due to collisions, and f2f is 
the term representing the color perturbation. The vectors ej are the velocity vectors 
along the links of the lattice. 

In this paper we use a triangular lattice (N = 6) for two-dimensional simulations 
with ej = (cos (27r(z — l)/6), sin (27r(z — l)/6)), and a body-centered-cubic lattice (N 
= 14) in three dimensions with G (±1, 0, 0), (0, ±1, 0), (0, 0, ±1) and (±1, ±1, ±1). 
For computational efficiency, we have used the single time relaxation model with 
the linear collision operator: 

^■ = -\{U-fl''^). (2) 
where r is the characteristic relaxation time, and f^'^'^^ is the local equilibrium distri- 
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bution given in two dimensions by 



and 



= ^ - ^>u^ (3) 



/^' = f + fe..u+f(e..u)^-^u^ (4) 



and in three dimensions by 



^(eq) ^ P _ P^2 



/i = g + gie* ■ u) + -(Gi • u) - gU , (6) 
for Gj along the lattice axes, and 

64 24^ ^ 16^ ^ 48 ^ ^ 

for ej along the links to the corners of the cube. In the above equations p(x, t) = 
X)j/i(x, t) and pu(x, t) = J2i fi{^,t)^i are the local density and momentum, respec- 
tively. 

The detailed forms of the coefficients in Equations (3-7) are determined by the 
conservation of mass and momentum, the constraints of Galilean invariance, and 



a velocity independent, isotropic pressure tensor. It can be shown that the 
macroscopics of (1-4) and (1,2,5,6,7) correspond to the incompressible Navier-Stokes 
equation in two and three dimensions, respectively: 

dtu + (u ■ V)u = — VP + z/V^u 

P 

V-u = 0, 

where P is the pressure, and u is the kinematic viscosity. In two dimensions u = 
(2r — l)/8, and in three dimensions u = (2r — l)/6. 



Defining the local order parameter as ipix,t) = J2iLo{fii'^,t) — t)), and the 
local color gradient G(x) by 

N 

G(x) = 5:e,{^(x + e,)}, (8) 

i=l 

we add the surface-tension inducing perturbation 

Q^i = A\G\cos2{di-dG) (9) 

to facilitate segregation and stabilize interfaces. Here 6i is the angle of lattice direction 
i, and 6g = arctan(Gj^/G^) is the angle of the local color gradient. The constant A 
sets the surface tension, a, through a ~ Arp. Note that G is perpendicular to 
red-blue interfaces and its magnitude large there, while in a homogeneous (color) 
region it approaches zero. In the recoloring step we then maximize the color flux 
H = J2f=iifi — fi)^i in the direction of the color gradient, G, by maximizing H ■ G. 
Recoloring conserves the individual color components and hence the total density of 



the fluid. It has been shown |jTO[ that Laplace's law holds for this model. 

For our simulations we carry out critical quenches with J2x V^(X) t) = 0. The largest 
systems we simulated were (2048)^ in two dimensions and (128)'^ in three dimensions. 
Although we have investigated the domain growth and scaling properties for a variety 
of lattice sizes and parameters, we report on the domain growth and dynamical scaling 
properties for only one set of parameters in both two and three dimensions. The 
results obtained with smaller lattices and different parameters (surface tension and 
initial fluctuations) are consistent with the data presented here. 

We initialize the lattice (in both two and three dimensions) with (^/'(x)) = and 
(u) = with small local fluctuations, where the angle brackets () signify an average 
over the lattice. The system then evolves according to the dynamics outlined above. 
In two dimensions the surface tension inducing parameter A = .01, and the average 
fluid density per site is p = 2.1, and in three dimensions A = .001 and p = 2.4. As the 



system evolves, single color domains form and grow while the total fluid undergoes 
a Navier-Stokes evolution. In Figure 1 we show a "snapshot" (for a system of size 
(256)^) of a two dimensional system at time t = 1000. 

One convenient way to characterize the growth kinetics during the segregation 
process is through the order parameter correlation function, G{r,t) = {ip{r)ilj{0)) — 
(ip)"^, averaged over shells of radius r. One can then deflne the domain size, R{t), as 
the flrst zero of the function G{r, t). The Fourier transform of G is then the structure 
factor S{k,t). As time evolves the structure factor becomes more sharply peaked, 
and its maximum moves to smaller values of the wavenumber, k. In a wide variety 
of phase segregating systems S{k, t) has been observed to follow the dynamic scaling 
relation at late times 

S{k,t) ^ R\t)F{x), (10) 

where x = kR{t), d is the spatial dimension and F{x) is the structure factor. 

In Figure 2(a) we plot R{t) for the two dimensional model. The data indicate 
that R{t) ~ t^^. This exponent is in excellent agreement with the generally accepted 
theoretical prediction of t^/^ [0 and the numerical simulations of Ferrel and Vails 
which flnd R{t) ~ t° where a ~ .65 and a ^ .69 for systems with and without 
thermal modes respectively. 

In Figure 3 we show the scaled structure factor, -F(x), for our two dimensional 
simulations. For x > 2 we flnd that F{x) ~ x~^'^ in reasonable agreement with 
Porod's law which predicts that F{x) ~ x"*-'^^^-* in this region. However, for x < 1, 
we observe that F{x) ~ a;^'^. This is in sharp contrast with the theoretical arguments 
of Furukawa @] and Yeung |15] which predict that F{x) ~ x^. Our flndings for x <^ 1 



appear to corroborate other recent numerical simulations of hydrodynamic spinodal 
decomposition in two dimensions. Shinozaki and Oono in their CDS model observe 
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that F{x) ~ x^. They conjecture that this might be a result of finite size effects 
coupled to fiuctuations of domain walls due to long range hydrodynamic interactions. 

In Figure 2(b) we plot R{t) for our three dimensional model. Here we find that the 
domain growth is approximately linear in time with R{t) ~ This is in excellent 
agreement with the theoretical predictions of Furukawa 0] and Siggia and with 
the (TDGL) numerical simulations of Ferrel and Vails and the CDS simulations of 
Koga and Kawasaki. It disagrees, however, with recent MD simulations of Ma et. 
al who observe a domain growth R{t) ~ ^-ssi.os^ ^j-^jg latter exponent may not 
refiect the growth in asymptotic time regime. 

In Figure 4. we plot F{x) for the three dimensional model. As in two dimensions 
we find good agreement with Porod's law with F{x) ~ x~^'^ for x > 3. Again, 
though we find marked disagreement with the small x predictions of F{x) which call 
for F{x) ~ x^. In particular we find that F{x) ~ x^'^. This result is consistent 
with recent light scattering data from spinodal decomposition of isobutyric acid and 



water |]T^ which also indicates a reduced exponent (~ 2), but differs from the results 
of Bates and Wiltzius ||I8| . 

We summarize by pointing out that the results of our simulations tend to be 
consistent with existing theories of domain growth in both two and three dimensions. 
There is also good agreement with Porod's law in both of these cases. However, we find 
a marked discrepancy with the theoretical predictions of Yeung [l^ and Furukawa 
for the small x behavior of the scaled structure factor. Moreover, there seems to be 
a consistent deviation from these theories among some recent numerical and some 



experimental work IjT^. In light of these results, we believe that the question of the 
small X behavior of the scaled structure factor is still open and that it is necessary to 
develop a more complete theory of dynamical scaling which includes hydrodynamic 
interactions. 
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The current model simulates fluids for deep (near zero temperature) quenches. 
As a result, this precludes an analysis of binary fluids near their critical points. To 
study critical dynamics of binary fluids requires some notion of temperature. An 
extension of the model described here might include a stochastic term which mimics 
the effects of thermal noise (Landau-Lifshitz fluctuating hydrodynamics) ||T^ or a 
kinetic temperature as in Reference with which one can control the segregation 
process by a local temperature. Without such a noise term it is unlikely that this 
model can simulate off critical quenches. 

The model in this paper is ideally suited for the simulation of hydrodynamic phase 



segregation in high Reynolds number flows . Moreover, the LB scheme can easily 
simulate spinodal decomposition in systems with complicated boundaries, stirring, 
and the effects of wettability. 
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FIGURE CAPTIONS 

Figure 1: Typical configuration for the two-dimensional model at time t = 1000. 
The system size is (256)^, and the other parameters are as in the text. In both (a) 
and (b) The black region represents sites with positive order parameter. 

Figure 2: Time dependence of average domain size for (a) two-dimensional and (b) 
three-dimensional systems. Both two and three dimensional data is averaged over 2 
independent runs. 
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Figure 3: Scaled structure function F{x) for 2 dimensions at various times: t = 
1000 (A); t = 2000 (□), and t = 3000 (O). The data is averaged over 2 independent 
runs. 

Figure 4: Scaled structure function F{x) for 3 dimensions at various times, t = 
1200 (A); t = 2400 (□), and t = 3200 (O). The data is averaged over 2 independent 
runs. 
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